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ABSTRACT 

In  reliability  growth  models,  systems  undergo  an  improvement  in  perfor¬ 
mance  during  prototype  testing,  as  design  changes  are  made,  and  operating 
procedures  and  environment  are  modified.  In  the  learning-curve  models, 
this  improvement  occurs  continuously  over  time,  and  there  is  great  in¬ 
terest  in  predicting  the  ultimate  performance  of  the  system,  using  only 
the  epochs  of  the  failures  which  occur  early  in  the  testing  program. 

This  paper  constructs  a  general  framework  in  which  to  analyze  this  prob¬ 
lem,  including  as  special  cases  many  different  model  variations  that 
have  previously  been  analyzed.  Numerical  trials  indicate  the  difficulty 
of  using  classical  procedures  to  estimate  ultimate  performance;  the 
maximum  likelihood  estimator  is  unstable  for  small  testing  intervals  with 
a  small  number  of  systems  on  test,  and  is  even  inconsistent  for  a  large 
number  of  systems.  Bayesian  procedures  are  recommended'  for  implementa¬ 
tion,  as  they  can  use  the  data  from  any  testing  protocol. 


A  GENERAL  FRAMEWORK  FOR  LEARNING  CURVE 
RELIABILITY  GROWTH  MODELS 

by 

William  S.  Jewell 


1.  INTRODUCTION 

When  complex  systems  are  first  produced  and  undergo  development  and 
testing  prior  to  actual  operational  use,  a  learning  curve  effect  usually 
takes  place,  in  which  the  initially  high  failure  rate  diminishes  over  time 
as  the  various  causes  of  breakdown  are  identified  and  corrected  through 
engineering  design  changes,  environmental  modification,  operating  pro¬ 
cedure  changes,  and  so  on.  This  "shakedown"  process  is  essential  to  the 
development  of  a  robust  system  that  can  be  routinely  produced,  installed, 
and  operated  with  high  reliability  under  long-duration,  adverse-environment 
circumstances.  Because  of  the  high  development  costs  of  such  systems,  it 
is  clearly  important  to  be  able  to  estimate  this  reliability  growth ,  and, 
especially,  to  estimate  the  ultimate  level  of  performance. 

A  large  variety  of  special-purpose  growth  models  have  already  been 
proposed;  [7]  and  [10]  are  convenient  summaries.  For  example,  one  group  of 
models  focuses  upon  the  defect-identification  process,  assuming  that  each 
system  failure  gives  the  engineer  a  chance  to  rectify  a  newly-discovered 
cause  of  failure.  Most  of  these  models  are  Markovian  in  nature;  see,  e.g., 
[2],  [3],  [6].  The  primary  drawback  of  these  models  seems  to  be  the  need 
to  specify  an  appropriate  defect  state  space,  and  to  estimate  the  related 
defect  removal  probabilities. 

A  second  method  of  modelling  reliability  growth  might  be  called  the 
structural  parameter  approach.  Learning  is  assumed  to  occur  only  in 


discrete,  well-defined  stages,  with  unknown  failure  parameters  associated 
with  each  stage;  only  the  number  of  failures  in  each  stage  is  used  to 
estimate  these  parameters.  If  growth  is  assumed  to  take  place  at  each 
stage,  then  this  leads  to  a  monotone  estimation  problem  in  which  the 
many  parameters  are,  by  hypothesis,  ordered  [2],  [21].  It  is  also  possible 
to  weaken  this  assumption  to  stochastically  ordered  parameters  [15],  [20]. 
The  most  serious  criticism  of  this  approach  to  reliability  growth  is 
that  it  gives  very  little  predictive  capability,  since  future  learning 
is  only  loosely  constrained. 

A  third  group  of  models  requires  strong  assumptions  about 
the  general  form  of  reliability  growth,  based  upon  past  experiences  with 
similar  processes;  then  one  estimates  a  relatively  few  free  parameters 
using  only  failure  epoch  data.  We  will  call  this  the  learning  curve 
approach,  since  it  is  closely  related  to  deterministic  models  of  produc¬ 
tion  learning,  developed  by  industrial  engineers  in  the  50's. 

A  large  literature  exists,  mostly  focused  around  special  growth  curve 
forms,  such  as  the  Duane  model  [7],  [10].  Many  of  the  recent  software 
reliability  models  [11],  [18],  [19]  are  also  of  this  type. 

The  objective  of  this  paper  is  to  present  a  general  framework  for 
constructing  and  estimating  the  parameters  of  learning-curve  models, 
in  particular,  two  important  variations  of  the  general  form,  the 
"as  operated,"  and  "as  produced"  models.  It  is  shown  that,  independent 
of  the  form  of  learning  deemed  appropriate,  the  type  of  model  selected 
as  well  as  the  testing  protocol  used,  the  parameter  estimation  problems 
are  mathematically  very  similar.  Preliminary  numerical  trials  indicate 
that  classical  point  estimators  of  the  parameters  may  be  misleading  as 
to  their  accuracy;  accordingly,  we  shall  outline  several  possible 


S' 


i- 


» 


h- 


u 


u 


**»  -#t-  ?-■  i". 


3 


Bayesian  approaches  to  parameter  estimation.  Detailed  results  from  our 
numerical  simulation  are  then  presented.  Our  main  conclusion  is  that 
most  learning-curve  testing  protocols  provide  little  of  the  parameter 
precision  needed  by  manufacturers  to  make  "early  returns"  predictions 
of  ultimate  reliability  using  traditional  estimation  methods. 


V  V  • 
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2.  FAILURE  MODEL  FORMULATION 

In  developing  a  general  model  of  reliability  growth,  we  must 
consider  carefully  the  assumptions  made  about  the  underlying  failure 
process.  First  of  all,  by  ignoring  any  specific  information  about  the 
type  of  failure  (except  to  assure  that  this  information  is  used  in  a 
positive,  predictable  way  by  the  engineers  to  modify  future  performance), 
we  restrict  ourselves  to  observing  point  processes .  Figure  1  shows  a 
single  system  on  test,  with  various  failure  events  ("glitches")  in¬ 
dicated  by  asterisks.  Our  objective  is  to  develop  a  stochastic  law 
governing  the  occurrence  of  these  events;  in  reliability  growth,  we  expect 
these  events  to  occur  less  frequently,  or  to  be  more  widely  spaced,  as 
time  evolves  from  the  test  origin.  (Section  7  discusses  results  with 
several  systems  on  test.)  Because  of  the  concerns  expressed  in  [14],  [16], 
a  failure-rate-oriented  model  appears  to  be  more  appropriate  than  an 
interval-oriented  one. 


-tt- 

dx 


x 

t 


lifetimes 
point  process 
failure  epochs 

local  time  (or  age) 
global  (or  clock)  time 


FIGURE  1.  Notation  with  Single  System. 


To  fix  notation,  let  . ..  be  the  successive  random  lifetines 

(intervals  between  failures),  and  t^,^.  •••  be  the  corresponding  failure 
event  epochs ;  clearly  tQ  =  0  ,  and  t^  ■  x^  +  x^  +  . . .  x^  (j  =  1,2,  ...)  . 


A  random  variable  is  indicated  with  a  tilde;  its  actual  value  in  some  reali¬ 
zation  is  the  same  variable  without  a  tilde. 
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Suppose  we  are  formulating  the  stochastic  law  governing  the  jth  lifetime. 

It  is  convenient  to  refer  to  a  small  interval  of  time,  dx  ,  which  is 
measured  either  in  local  time  (or  age) ,  x  ,  since  the  actual  last  fail¬ 
ure  epoch,  tj_^  ,  or  in  global  (or  clock)  time,  t  =  t^  ^  +  x  ,  from 
*  th 

the  test  origin.  A  failure-rate  law  for  the  j  interval  would  then  be 
completely  determined  if  we  could  specify  the  probability  of  the  next  failure 
epoch  occurring  in  dx  as  a  function  of  the  entire  past  history  of  the  point 
process,  i.e.,  as  a  function  of,  say,  (x^.x^,  •••,  x^_^;x)  ,  or,  of 

(t^tj,  ...,  t^_1;t)  . 

Our  first  assumption  will  be  that  the  starting  epoch  of  an  interval 
is  a  surrogate  for  all  previous  failure  history  of  the  system;  i.e.,  for 
the  interval,  x^  and  t^  are  statistically  dependent  only  on  t^  ^  = 
t^_^  ,  and  not  on  the  other  failure  epochs  (and  not  upon  the  index  j). 
Physically,  this  means  that  the  global  time  t^_^  describes  sufficiently 
all  of  the  engineering  changes,  the  characteristics  of  the  replacement 
items,  the  procedural  changes  that  have  been  incorporated  into  the  system 
to  date,  etc.  Mathematically,  we  are  specifying  a  failure  rate  function, 

irk 

f  ,  such  that: 


(2.1) 


Pr  {xj  e  (x,x+dx]  |  x.  >  x  ;  x^Xj,  ...  x^^it) 
=  f (x; t  |  tj_x  «  x1  +  x2  +  ...  +  x^_1)dx 


for  (j  =  1,2,  ...)  ,  (0  <_  x  <_  Xj)  ,  and  (t^  ^  <  t  <_  t^)  .  For  simplicity 
in  the  sequel,  we  will  move  freely  between  t  and  x  *  t  -  t^  ^  ,  whenever 
the  interval  index  j  is  clear  from  the  context. 


Time  and  age  can  also  be  measured  in  production  time  units,  or  in  actual 
system  operational  time.  This  raises  certain  measurement  difficulties, 
but  no  conceptual  ones. 

t 

We  emphasize  the  case  where  the  random  variables  are  continuous,  but 
discrete  cases  follow  directly,  with  only  minor  modifications. 
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The  second  assumption  we  shall  make  is  that  the  age-dependent  and 
time-dependent  failure  meshanisms  of  the  system  are  statistically  in¬ 
dependent ,  e.g.,  that  the  failure  rate  function  is  the  sum  of  two  terms, 
so  that  for  interval  j  , 

(2.2)  f(x;t  |  t  x)  =  <j>g(t  |  t^  j)  +  mh(x)  . 

Here  <|>  and  to  are  two  parameters  to  be  estimated,  and  g  and  h  are 
two  prototypical  growth  forms,  assumed  known  from  previous  experience, 
h  is  the  usual  age-dependent  hazard  function,  while  g  is  the 
learning-curve  itself,  which  is  decreasing  in  some  given  manner  as  t 
increases. 

This  assumption  is  reasonable  if  we  think  of  a  single  unreliable 
item  being  replaced,  with  the  "wear-out"  hazard  h  being  a  property 
of  this  item,  while  g  reflects  independent  improvement  of  the  rest 
of  the  system  or  the  operational  environment.  If  the  system  consists 
of  many  unreliable  terms  with  exponential  lifetimes  (h(x)  =  1)  ,  then 
this  model  also  applies,  if  the  reliability  growth  provides  independent 
environmental  improvement,  or  if  it  rectifies  or  removes  certain  of  the 
unreliable  elements  [5],  [18].  Further  development  or  validation  of 
(2.2)  under  more  complex  system  assumptions  is  needed. 
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3.  ALTERNATE  FORMULATIONS;  "AS  OPERATED"  AND  "AS  PRODUCED"  GROWTH  MODELS 

V 

The  choice  of  the  learning  curve  function,  g  ,  depends  upon  further 

*  assumptions  about  the  nature  of  the  improvements  made  upon  the  system 
over  time.  We  shall  consider  two  special  cases  which  we  believe  bracket 
the  possibilities  to  be  found  in  practice. 

Type  I  ("as  operated")  models  refer  co  those  situations  in  which 
the  reliability  growth  is  either  due  to  improvements  in  environmental 
conditions  (cooling,  shock  insulation,  power  supplies,  support  equipment 
etc.),  continuing  changes  in  operational  procedure;,  or  when  any  design 
improvements  on  current  production  models  can  be  retrofitted  on  existing 
systems  in  the  field  or  on  test.  In  this  case,  we  assume  that  a  one- 

5 

dimensional  prototypical  learning  curve  g(t)  is  known,  and  that  in 

.  (2.2): 

(3.11)  g(t  |  tj_x)  =  g(t)  ,  (t^  <  t  <_  tj  =  tj_x  +  x^)  . 

Most  models  in  the  literature  are  of  this  type. 

Type  II  ("as  produced")  models  will  refer  to  the  other  extreme  in 
which  reliability  improvements  are  built  into  each  system  only  during  the 
production-installation  process,  and  thus  do  not  affect  equipment  already 
in  the  field  or  on  test.  As  above,  there  is  an  underlying  learning  curve 

g(t)  ,  but  for  the  interval,  only  the  value  of  the  function  at  the 

starting  epoch  is  of  interest,  and 

*  (3. Ill)  g(t  |  t  x)  =  g(tJ_1)  (tJ_1  <  t  <  tj) 

«• 

is  a  constant  over  this  lifetime  interval.  Models  of  this  type  were 
introduced  in  [4],  [5],  [18].  Note  that  a  Duane  growth  curve  (see  below) 
has  some  technical  difficulty  in  the  first  interval  of  model  II. 
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One  can  easily  imagine  combinations  of  these  two  types  of  models, 
as  when  retrofitting  engineering  changes  gives  less  improvement  when 
applied  to  used  items.  However,  Models  I  and  II  are  clearly  good 
approximations  to  each  other  when  g  is  a  slowly  varying  function, 
as  compared  to  the  random  lifetimes,  so  that,  in  many  situations,  we 
may  not  expect  these  models  to  give  dramatically  different  results. 

Figure  2  shows  how  the  actual  failure  curves  for  the  two  different  models 
would  be  constructed. 


r<  t*. 


FIGURE  2.  Components  of  failure  rates,  uih  and  <f)g  ,  and 
resulting  total  failure  rates,  f  ,  of  Models  I 
and  II  for  a  particular  realization. 


4.  CHOICE  OF  GROWTH  FUNCTIONS 


There  remains  the  problem  of  specifying  the  form  of  the  learning 
curve,  g  .  Much  of  the  current  literature  specifies  Crow's  adaptation 
of  a  deterministic  learning  curve  of  Duane  [1],  [7],  [8],  [9],  [10],  [17], 
in  which  f(x,t  |  t^  =  g(t)  *  ktT  ,  with  the  internal  parameter 
Y  <  1  for  reliability  growth.  The  resulting  point  process  is  referred 
to  as  a  Weibull  process.  In  our  view,  this  model  is  incomplete  because 
there  is  no  residual  failure  effect  (w  «  0)  at  t  00  .  Also,  this  curve 
"blows  up"  for  t  =  0  ,  which  makes  the  normalizations  described  below 
in  Sections  5  and  9  somewhat  awkward.  We  would  prefer  to  use  a  more  tract- 
able  form,  such  as  g(t)  =  e  ,  for  example.  However,  much  practical 
work  remains  in  the  selection  and  justification  of  a  particular  learning 
curve,  since  it  seems  extremely  difficult  to  identify  the  correct  form 
of  g  with  the  usual  amounts  of  data  available. 

The  "hidden"  shape  parameter,  y  ,  in  any  g  is  also  a  problem. 
Normally,  one  would  like  to  leave  this  and  similar  parameters  to  be  estimated 
by  the  procedure  described  below.  However,  these  parameters  are  quite 
non-linear  in  their  effect  on  the  likelihood,  and  in  the  few  numerical 
trials  we  have  made,  are  extremely  sensitive  to  small  variations  in  the 
data. 

It  is  also  clear  that  nothing  prohibits  us  for  specifying  a  g  which 
changes  value  only  at  certain  pre-defined  points  in  time,  rather  than  con¬ 
tinuously.  But  perhaps  such  situations  are  better  handled  through  the 
discrete  structural  parameter  models  described  in  Section  1. 


5.  DISTRIBUTIONS  AND  DATA  LIKELIHOODS 


The  random  lifetimes,  x^,X2>  ...  , 
dependent  and  not  identically  distributed 
and  a))  the  conditional  distribution  of 


generated  by  (2.2)  are,  in  general, 
.  However,  given  (and  4> 

Xj  is  found  from  the  relation 


(5.1)  f(t,x  |  tj^)  -  -  in  PC(x  |  tj^j^.uj)  , 

c 

where  P  is  the  complementary  distribution  function  defined  by 


(5.2)  PC(x  |  -  Pr  {Xj  >  x  |  tj_1;<j»,u)  , 

for  (j  *  1,2,  ...)  (x  >. 0)  and  (t  ■  t^_^  +  x)  .  The  density  of  x^  is 

c 

just  the  negative  derivative  of  P 

Let  G(t)  and  H(x)  be  the  integrals  of  the  growth  and  age  hazard 
components,  measured  from  the  test  origin  and  the  current  interval  origin, 
respectively: 


(5.3) 


G(t) 


g(u)du 


;  H(x) 


h(w)dw  . 


(As  described  above,  this  raises  minor  technical  difficulties  in  the  Duane 
model  in  which  G(0)  is  undefined.)  Then  it  is  straightforward  that  for 
model  type  I: 

(5.41)  PC(x  j  tj_1;<t>,uj)  -  exp  (-<t>[G(t)  -  G(t  -  oiH(x)}  , 

whereas  for  type  III 


(5.4II)  PC(x  |  tj_^ ;4> ,co)  -  exp  {-<(>xg (t^_1)  -  uH(x)}  , 

remembering  always  that  t  ■  t.,  +  x  in  the  current  interval  j  . 


Having  developed  the  law  for  each  interval  of  our  point  process,  we 

can  now  consider  specific  testing  protocols,  calculate  appropriate  data 

likelihoods,  and  proceed  with  the  estimation  of  <J>  ,  u  ,  and  any  hidden 

parameters  in  g  .  Remember  that  we  have  assumed  g  and  h  are  known. 

Suppose  there  is  one  system  on  test,  as  in  Figure  1,  and  the  test 

is  stopped  at  some  predetermined  global  time  t  ■  T  ("time-truncated"). 

During  (0,T]  ,  there  will  have  been  a  random  number  n(T)  =  sup  {n  |  t  < 

of  complete  lifetimes,  x^Xj,  ...,  x~^  »  together  with  the  (h(T)+l)St 

"in  progress"  lifetime  with  current  age  x  *  T  -  t~.„N  .  The  appropriate 

n(T; 

data,  V  ,  from  our  point  process  will  then  be: 

(5.5)  V  =  (n;T;x. ,x. ,  ...,  x  }  . 

i  /  n 


The  likelihood  is: 

n 

(5.6)  L(<t>,u  |  V)  =  n  [*g  (D)  +  coh.(P)]  exp  {-*G(P)  -  ujH(P)}  , 

j=l  2  3 

where  for  both  models  I  and  II: 

n 

(5.7)  h  (V)  =  h(x.)  ;  H(V)  =  J  H(x.)  +  H(T  -  t  )  . 

j  j  j=l  j  n 

H  is  called  the  "total-Q-on-test"  in  [12]. 

The  learning  curve  data  components  depend  on  the  model  type;  for 
growth  "as  operated": 

(5.81)  gj (V)  =  g(t  )  ;  G(V)  -  G(T)  ; 


while  for  growth  "as  produced": 


(5.8II)  gj(^)  =  gCt^)  ;  G(V)  =  l  Xj8<tj-l)  +  _tn^8^tn^  ‘ 


It  is  clear  that  these  components  are  quite  close  for  slowly  varying  g  . 

If  observations  are  stopped  at  the  n*"*1  failure  epoch  ("failure  truncated") 
the  above  formula  still  apply,  with  T  =  t 
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The  classical  maximum-likelihood  estimates  <t>  and  id  follow  directly 
from  (5.6);  we  get  an  implicit  solution  through  the  two  equations: 


(6.1) 


n  h. 

I  — 1 — 

j-l  $gj  +  whj 


n  8 . 

H  ;  l  — - 

j-l  <*)gj  +  oihj 


If  g  also  contains  an  internal  parameter  y  ,  then  we  have  at  y  =  y 


(6.2) 


3g, (y) 


j-l  (Jigj  +  whj 


3G(y) 

3y 


in  an  obvious  extension  of  notation.  For  computational  convenience,  we 
record  also  the  second  derivatives  of  i.n-likelihood  at  the  maximum  likeli- 

AAA 

hood  point  (<j>,w,y)  : 


(6.3) 


a2  i  n  (g,)k(h,)2"k 

k  ^ lk  -  -  I  -T1 - 5 - 2  ;  (k  =  °‘1‘2) 

a*  3oi  j-l  (4>gj  +  uh^T 


(6.4) 


1  3^  £n  L 


j-l  (4>gj  +  uihj)' 


I  3  L 
$  3^3y 


3 2  fcn  L 


2  .  .  2 
3  g ,  /3gA 

n  (*8j  +  “V  —2-  -  *{-£) 


3Y  j-l 


(4>gj  +  whj) 


Because  these  derivatives  are  evaluated  at  [ 4> , to , y  1  »  they  are  also  equal 
to  L  1  times  the  corresponding  derivatives  of  L  . 


OTHER  TESTING  PROTOCOLS 


The  results  of  Section  6  can  be  easily  adapted  to  other  testing 
arrangements.  For  example,  suppose  that  S  systems  were  put  on  test 
at  the  same  global  time  origin,  and  operated  for  the  same  total  interval 
T  .  Then,  system  i  would  report  n^  failures,  with  completed  lifetimes 
Xil’Xi2’  Xin  ’  anc*  t*ie  tota^  data  set  corresponding  to  (5.5)  would 


(7.1) 


V  i=l  jni;T;Xil’Xi2’  **■’  Xinij 


with  (5.7)  replaced  by: 


(7.2)  hi;j(0)  *  h(x  )  ;  H±(V) 


X  H(V  +  H(T  - 


g^ (V)  and  G^(V)  are  defined  analogously  to  (5.81)  or  (5.8II);  then 
the  total  L(V)  is  the  product  of  terms  like  (5.6)  for  each  system 
i  =  1,2,  ...»  S  . 

It  follows  that  the  maximum  likelihood  estimator  formulae  (6.1)  and 

(6.2)  become: 


S  ni 


(7-3)  l  l 


i=l  j=l  <J>g±j  +  wh 


l  h  ;  l  l 

i=i  1  i=i  j= 


s  “i  g, 


i=l  j=l 


l  h  s 

i=l  1 


(7.4) 


S  "i 

I  l  ■:  ay  , 

i=*l  1=1  <(>g.  .  +  oih ,  , 


S  3Gi( Y) 
i=l 


A  variety  of  other  testing  protocols  can  easily  be  incorporated  into 
similar  formulae,  including  different  test  origins  and  terminations, 
missing  intervals,  etc.;  the  stopping  time(s)  may  also  depend  upon  the 
data  outcome  (see,  e.g.,  [12],  [13]). 


8.  CONSTANT  AGE  HAZARD 


It  is  interesting  to  note  that  only  the  ratios  g^/h^  entered  into 
the  computations  (6.1)  -  (6.5);  in  effect,  one  could  switch  to  an  operational 
time  and  analyze  an  equivalent  model  in  which  h(x)  =  1  for  all  x  . 

This  constant  age-hazard  is  important  in  ordinary  reliability  applica¬ 
tions,  where  it  corresponds  to  exponential  lifetimes  and  the  stationary 
Poisson  failure  process.  It  is  of  interest  to  see  what  happens  in  the 
corresponding  reliability  growth  model,  with  h^  =  1  for  all  i  ,  and 
H  =  T  . 

Models  of  type  I  then  give  rise  to  time-varying  (inhomogeneous) 

Poisson  processes,  with  intensity  parameter 


(8.1) 


Mt)  =  <J>g(t)  +  a)  . 


The  distribution  of  the  number  of  failures  in  (0,T]  is  Poisson: 


(8.2)  Pr  (h(T)  =  n  |  <f>,w}  =  [cJ>G(T)  +  uT]ne_[l)>G(T)+alTl/n! 


from  which  we  can  get  the  law  of  n  failure  epoch  through  the  equivalence 

[t  <  T]  [n(T)  >  n]  .  The  mean  and  variance  of  n(T)  are,  of  course, 

n  —  — 

<J>G(T)  +  o)T  .  Moments  of  the  lifetimes  {Xj}  ,  on  the  other  hand,  have 
no  simple  form  in  the  general  case;  however,  the  concept  of  looal  mean 
life  at  t  ,  introduced  in  [4],  and  defined  as: 


(8.3) 


y(t)  =  1/A(t)  =  [<J>g(t)  +  u]“A  , 


has  the  interpretation  of  mean  duration  until  the  next  failure  if  no 
further  change  occurs  in  g(t)  (or  if  it  is  slowly  varying).  Of  course, 
if  g(t)  0  as  t  00  ,  successive  lifetimes  approach  independent. 


exponentially-distributed  random  variables  with  mean  life  y(“>)  =  w 


.<  ik 


Models  of  type  II  with  h(x)  =  1  ,  on  the  other  hand,  give  rise  to 

independent,  exponentially-distributed  random  lifetimes  for  every 
t  h 

interval;  the  j  lifetime  has  mean: 

E{x^  |  tj_1; 4>,to}  =  vKt^)  =  [<frg(tj_1)  +  u>]  1  . 

Note  that  this  point  process  is  not,  however,  a  Poisson  process,  since 
the  intensity  parameter  jumps  to  its  new  value  only  at  a  failure  epoch, 
and  not  continuously  over  time.  Hence,  there  is  no  simple  form  for  the 
distribution  of  n(t)  . 
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9.  BAYESIAN  ESTIMATION 

Because  most  product  development  programs  Involve  systems  which  are 
evolutions  of  previous  ones,  it  makes  sense  to  quantify  this  engineering 
experience  and  "know-how"  by  using  Bayesian  estimation  procedures.  There 
appear  to  be  three  major  ways  in  which  this  approach  might  be  applied 
in  reliability  growth. 


A.  Full  Bayesian  Estimation 

In  a  full  Bayesian  analysis,  the  likelihood  (5.6)  (or  an  extended 
version  developed  in  Section  7)  is  used  together  with  a  joint  prior  density 
p (4> , o>)  in  the  usual  Bayes'  formula  to  find  the  posterior  joint  density : 


(9.1) 


p(<f>,u>  |  V ) 


L(<t>,o)  1  P)p(<fr,m) 
normalization  (P) 


There  is  probably  little  to  be  gained  by  attempting  an  analytic  solution 
to  this  equation  through  the  use  of  natural  conjugate  priors.  Even  if 
p(<j>,u))  factored  into  independent  gamma  densities  on  <f>  and  o>  ,  it  can 
be  seen  from  (5.8)  that  the  posterior-to-date  density  would  have  n 
terms,  corresponding  to  different  priors  ^  (j  =  0,1 . n)  . 

If,  in  addition  to  $  and  w  ,  a  hidden  parameter  in  g  or  h 
must  be  estimated,  this  leads  to  a  very  non-exponential-type  likelihood, 
for  which  no  natural  conjugate  prior  exists. 

Thus,  in  general,  it  seems  that  useful  full  Bayesian  solutions  of 
(9.1)  would  require  2-  or  3-dimensional  numerical  integration. 

B.  Multinormal  Approximation 

If  the  prior  and  the  likelihood  are  appropriately  bell-shaped 
functions  of  w  and  <f>  ,  then  one  can  make  multinormal  approximations 


I  1  Wf  ; 
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of  the  densities,  and  then  carry  out  an  exact  calculation  using  (9.1). 

The  parameters  for  the  likelihood  would  have  to  be  calculated  numerically 
using  the  formulas  in  Section  6;  the  multinormal  prior  would  require 
eliciting  the  prior  means  and  covariances  of  oj  and  $  from  engineering 
personnel,  or  examining  failure  rate  records  of  related  development 
programs.  At  the  very  least,  this  requires  knowing  how  much  data  and 
prior  confidence  is  needed  to  obtain  unimodal  densities  of  the  correct 
form. 


C.  Conditional  Likelihood 


The  estimation  problem  is  greatly  simplified  if  we  can  assume  that 
only  one  parameter  is  imprecisely  known.  A  possible  situation  in  which 
this  might  occur  is  where  the  shape  of  the  learning  curve  and  the  initial 
failure  rate  f(0,0  |  0)  =  a  are  well-known  as  testing  begins.  Assuming 
that  the  prototypical  failure  functions  are  normalized  so  that  h(0)  = 
g(0  |  0)  =  1  ,  we  set  a  =  (j>  +  w  ,  and  use  the  one-dimensional  likelihood 
and  prior  in: 


(9.2) 


p(u)  |  V,a ) 


L(u>  1  P,ct)p(o)  1  a) 
normalization 


S  °i  (  S  S 

(9.3)  L(io  |  V,a )  =  n  n  (ag  +  w(h  .  -g  .))  exp  !-a  £  G  -w  l  (H  -G  ) 

i=l  j=l  13  13  3  (  i=l  i=l  1  1 

Thus,  attention  is  shifted  to  estimation  of  the  ultimate  failure  rate, 
w  ,  which  is  usually  the  performance  parameter  of  greatest  interest. 

It  shall  be  remembered  though,  that  a  (and  y)  are  assumed  known. 

For  instance.  Figure  3  shows  a  typical  prior  density  p(io  |  a)  ,  with 
experience  indicating  a  most  likely  improvement  of  about  40%  from  the 
initial  failure  rate  a  ,  but  with  allowances  being  made  for  the  possi¬ 
bility  of  other  outcomes,  even  oj  >  a  (reliability  loss).  This  prior 


would  then  be  multiplied  by  L(u>  |  P,a)  in  the  usual  way  to  give  the 
(unnormalized)  posterior  p(u>  |  V,a)  .  The  point  estimate  of  u>  would 
then  depend  upon  the  loss  function  selected. 


FIGURE  3.  Illustration  of  Bayes'  Law  for  Conditional 
Likelihood  Approach. 
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10.  NUMERICAL  TRIALS  -  EXPONENTIAL  GROWTH  FUNCTION 
A.  Random  Variates 

To  illustrate  the  numerical  difficulties  associated  with  reliability 
growth  models,  a  series  of  random  variates  were  generated  using  constant 
hazard  rate,  an  exponential  growth  function  (assuming  Model  I  is  operating), 
and  a  generous  set  of  parameters  in  which  a  significant  decrease  in  failure 
rate  occurs  after  about  15-20  failures  have  occurred  in  each  realization. 

More  specifically,  pseudo  random  lifetimes  {*y  ;  i  -  1,2,  ....  32  ; 
j  =  1,2,  ...,  80}  were  generated  using: 

h(x)  =  1  a  =  1.0  ;  oj  =  0.2  (0  =  0.8) 

g(t)  =  e  Yt  y  =  0.05  . 

This  gives  a  time-varying  Poisson  process  in  which  the  failure  rate  decreases 
from  an  initial  rate  of  1.0  towards  an  ultimate  rate  of  0.2,  with  a  time 
constant  of  y  1  =  20  time  units. 

Figure  4  shows  the  resulting  32  point  process  realizations  over  the 
range  0  _<  T  <_  40  ,  corresponding  to  zero  to  two  time  constants,  which  is 
probably  the  range  of  practical  interest,  as  only  13.5%  of  the  growth  is 
unobserved  at  T  =  40  .  (Because  of  quantization  in  the  printing  of 
Figure  4,  failures  that  are  closer  together  than  0.5  time  unit  show  as 
only  one  event;  for  example,  In  realization  #1,  the  failures  near  T  =  2.5  , 
5.0  ,  and  14.5  are  actually  double  events.)  Even  with  the  high  variability 
of  the  Poisson,  the  strong  reduction  in  failure  rate  with  time  can  be 
visually  inferred. 

The  average  number  of  events  in  one  realization  follows  the  law: 

M(t)  =  E{n(t)}  =  Q.2t  +  16(1  -  e  ’^C)  . 
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Figure  5  is  intended  to  show  how  this  form  is  better  and  better  approxi¬ 
mated  by  £  n^(T)/S  ,  as  the  number  of  realizations,  S  ,  is  increased. 
First  of  all,  for  S  =  1  (realization  //l) ,  this  sample  function  is 
just  the  usual  unit-step  counting  function,  with  jumps  at  the  t^  ;  it 
starts  out  rather  well,  but  then  lags  behind  the  mean,  beginning  with 
the  rather  long  8th  lifetime.  When  the  second  realization  is  included, 
making  S  =  2  ,  the  jumps  in  the  sample  function  are  of  height  1/2  ,  and 
occur  at  failure  epochs  {t^  ;  t^j^  systems  on  test.  The  result¬ 

ing  curve  is  somewhat  smoother,  but  of  course  still  contains  the  bias 
of  the  results  from  system  #1.  Keeping  in  mind  the  jagged  nature  of  the 
actual  sample  functions,  the  remaining  curves  show  how  closely  the  sample 
function  approaches  M(t)  for  S  =  4,8,16  ,  and  32  ,  for  T  =  0(5)40  . 

The  S  “32  curve  deviates  from  M(t)  by  less  than  0.26  units. 


B.  Maximum-Likelihood  Estimates 

Maximum-likelihood  estimates  of  to  were  then  evaluated  using  the 
(9.3)  version  of  the  equations  of  (7.3),  that  is: 


S  “iVi'  h  -  g  S 

(10.1)  l  l  - ^ ^ -  =  l  (H  -  G  )  . 

i=l  j=l  agij  +  w(hij  -  g^)  i-1 

In  our  case,  of  course,  h^  =  1  ,  g^  =  exp  (-yt^)  ,  and  the  RHS 
is  just  S[T  -  G(T)j  .  These  values  were  obtained  through  an  interactive 
numerical  search  over  successively  finer  values  of  to  ;  when  the  maximiz- 

A 

ing  <o  was  obtained  to  within  0.01,  parabolic  interpolation  using  three 
neighboring  values  of  u  was  used  to  estimate  the  final  to  . 


FIGURE  5.  Sample  functions  for  M(T)  from  simulations  with  S  =  1,2,4,8,16,32  systems.  Only 

selected  values  of  T  shown  for  larger  values  of  S.  (Straight  lines  added  for  continuity. 
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For  small  testing  intervals,  T  ,  the  formal  solution  of  (10.1)  is 
often  negative;  these  values  have  been  retained  in  Figures  6,  7,  and  8, 
even  though  the  likelihood  is  then  monotone  decreasing  over  (0,“)  , 
and  theoretically  we  should  set  =  0  .  However,  negative  values  of 
a)  give  us  some  idea  of  the  resulting  steepness  of  the  likelihood  near 
zero,  as  well  as  the  algebraic  instability  of  (10.1)  for  small  T  . 

A 

However,  the  negative  values  of  oi  are  not  especially  exact,  as  they 
were  determined  by  parabolic  extrapolation  from  m  =  0,  .01,  and  0.2  . 

Figure  6  illustrates  the  behaviour  of  m(T)  as  a  function  of  T 
for  S  =  1  (realization  #1) .  Notice  the  typical  behaviour  of  point  process 
estimators;  jumps  in  value  occur  at  the  failure  epochs,  followed  by  gradual 
decreases  (with  decreasing  slopes  as  the  number  of  samples  increases)  as 
the  testing  interval  increases  without  an  intervening  failure.  For  T 
from  0  to  t^  ,  a)  is  undefined;  from  t^  to  about  T  =  7  ,  it  shows 
reliability  decay  (ai  >  1)  ,  and  from  T  =  7.5  to  about  27.5,  there  are 
intervals  of  negative  values.  Even  though  we  know  w(T)  -*■  w  almost  surely 
as  T  00  ,  realization  it  1  is  still  estimating  a  value  of  w  less  than 
0.01  at  T  «  50  (yT  -  2.5)  . 

Different  results  will  of  course  be  obtained  from  different  realizations. 
Figure  7  compares  the  results  obtained  using  the  first  four  realizations. 

Only  selected  values  of  T  =  5(5)50  are  evaluated;  the  actual  estimator 
behaviour  between  these  points  looks  like  Figure  6.  Here  we  see  the 
large  variability  in  behaviour  between  samples,  especially  for  yT  <  1.0. 
Samples  it 2,  it 3,  it 4  are  closer  to  the  true  w  at  yT  ■  2.5  than  it  1,  but 
there  is  still  some  variability  left  at  this  point,  where  the  numbers  of 
samples  in  each  realization  are  19,  19,  28,  and  26,  respectively. 


FIGURE  6.  Maximum-likelihood  estimator,  w  ,  for  S  =  1  ,  using  realization  //I 
versus  testing  interval,  T  .  (See  text  regarding  negative  values.) 


FIGURE  7.  Maximum-likelihood  estimator,  w  ,  for  S  =  1  ,  for  realizations  it 
for  selected  intervals,  T  .  (Straight  lines  added  for  continuity 


Clearly,  the  estimator  will  be  more  stable  if  we  increase  S  ,  the 
number  of  systems  on  test.  Figure  8  shows  u>(T)  versus  T  =  5(5)50 

for  S  =  2,4,8,16  ,  and  32  .  S  =  2  still  has  a  large  region  of  negative 

results  (because  to  of  both  realizations  Hi  and  #2  in  Figure  7  are 
negative),  but  S  =*  4  and  higher  are  relatively  stable  after  T  =  20 
(yT  =  1.0)  .  Remember  that  the  complete  behaviour  between  selected  T 
should  look  more  like  Figure  6,  but  that  with  increasing  S  ,  the  jumps 
will  be  smaller,  and  will  occur  more  often.  Note  also  that  the  vertical 

scale  of  Figure  8  is  finer  than  that  of  Figures  6  and  7,  so  that  the 

relative  stability  is  better  than  it  appears.  The  numbers  of  samples  at 
T  *  50  are  38,  92,  195,  395  ,  and  787  for  S  =  2,4,8,16  ,  and  32  , 
respectively. 

C.  Data  Information  Estimates 

For  the  cases  in  which  an  internal  maximum  occurs  in  the  likelihood, 
it  is  convenient  to  have  a  measure  of  how  sharp  this  maximum  is.  A  con¬ 
venient  choice  is  the  (Fisher)  Information,  defined  as: 


If  I  as  a  function  of  to  is  approximated  by  a  normal  curve  about  to  , 
then  I  is  approximately  the  precision  (inverse  variance)  of  the  approxi¬ 
mation.  Estimates  of  the  information,  I  ,  were  obtained  during  the 
course  of  the  computations  of  to  ;  if  the  formal  solution  of  (10.1)  gave 

A 

negative  values,  I  was  computed  at  these  maxima,  but,  of  course,  the 
information  then  has  no  useful  physical  interpretation. 

Figure  9  shows  the  values  of  I  obtained  at  selected  values  of 
T  -  5(5)50  for  different  values  of  S  .  Points  with  negative  u)  are 


FIGURE  8.  Maximum-likelihood  estimator,  w  ,  for  S  =  2,4,8  ,  and  16  systems, 

for  selected  testing  intervals,  T  .  (Straight  lines  added  for  continuity. 
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indicated  differently  than  those  with  positive  u  .  For  S  =  1  ,  the 
first  four  realizations  lead  to  quite  different  results  (and,  of  course, 
the  actual  curve  versus  T  is  discontinuous).  However,  as  S  increases, 
we  may  expect  the  sampling  variability  to  decrease,  and,  in  fact,  the 
curves  for  S  greater  than  four  are  quite  smooth,  and  more  stable  than 
those  for  u>  .  Note  particularly  that  the  curves  appear  almost  to  be  the 
same  curve,  displaced  by  the  same  amount  as  S  is  doubled. 

D.  Data  Likelihoods 

To  give  some  idea  what  the  different  precisions  mean  for  our  example, 
Figure  10  shows  the  likelihood  as  a  function  of  u  ,  p(P  |  m)  *  L(u  |  V)  , 
for  T  =  20  (yT  =1.0)  ,  and  various  values  of  S  .  For  S  =  1  and  2  , 
the  formal  solution  for  o>  gives  negative  values,  and  hence  the  likelihood 
is  monotone  decreasing  over  (0,°°)  ;  with  different  realizations  chosen, 

different  results  might  have  been  obtained.  For  S  =  4  and  above, 

/\  A 
internal  u>  are  obtained,  and  the  information  7  increases  from  about  24 

to  196,  as  the  likelihoods  become  sharper.  For  S  =  2  and  higher,  there 

is  practically  no  likelihood  that  w  is  greater  than  1.0  (reliability 

decay);  for  S  =  8  and  higher,  there  is  practically  no  chance  that  u> 

is  larger  than  0.7.  However,  in  all  cases,  the  instability  of  the  point 

estimate  w  is  clearly  shown,  and  there  is  still  a  rather  large  confidence 

region,  even  for  S  ■  32  (455  samples). 
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11.  ANALYTIC  RESULTS  ON  ESTIMATORS 


Although  the  numerical  results  of  the  last  section  were  enlightening, 
they  were  obtained  for  particular  parameter  values,  and  a  particular 
choice  of  learning  cut  te  form.  It  is  desirable  to  also  ask  what  addi¬ 
tional  results  can  be  obtained  analytically,  apart  from  the  well-known 
result  that  u)(T)  -*■  a)  almost  surely,  as  T  -+■  00  .  Unfortunately,  the 
complicated  form  of  (10.1)  permits  only  a  few  explicit  results. 

A.  Initial  Behaviour  for  S  =  1 

The  poor  initial  behaviour  of  the  single-system  estimators  can  be 
explained  by  noting  first  that,  for  T  <  t^  ,  no  estimators  exist.  For 
t^  T  <  t^  »  (10.1)  has  only  one  term,  so  that: 

(11.1)  w(T)  =  (Hx  -  G1)‘1  -  ag11(h11  -  gu)-1  , 
and 

(11.2)  I(T)  -  (Hx  -  Gx)2  . 

In  particular,  for  the  forms  assumed  in  the  example,  and  assuming  yt-Q 
is  small,  we  find: 

(11.1*)  m(T)  K  (2/yT2)  -  a(l/Ytu)  , 

(11.2*)  I (t)  «  y2t4/4  . 

This  means  that  the  initial  MLE  estimator,  uKt^)  »  has  the  following 
sensitivity  to  the  first  lifetime: 
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if 

“'a  * 2  • 

then 

“(tU)  <  0 

(11.3)  if 

2(1  -  (2y/ a) )  <  at^  <  2 

,  then 

o  <  m(tu) 

and  if 

atii  <  2^1  “  (2y/°0)  » 

then 

m(t^)  >  a 

(assuming  that  y  is  small  compared  to  a) .  Since  only  the  middle  estimate 
corresponds  to  reliability  growth,  this  means  that  the  initial  estimates  of 
co  are  likely  to  be  pathological  ones,  since  the  range  of  allowed  t^  is 
small  (1.9  <  t^  <2.0  in  our  example). 

The  estimator  of  the  precision,  on  the  other  hand,  is  a  smooth  function 
of  T  ,  and  indicates  how  quickly  precision  builds  up  after  T  >  t^  ,  un¬ 
til  the  next  discontinuity  at  T  *  t^  • 

B.  Initial  Behavior  for  Arbitrary  S 

The  above  discussion  can  be  extended  to  an  arbitrary  number  of  systems, 
provided  we  limit  our  discussion  to  values  of  T  between  the  first  failure 
to  occur, 

min  (tj^jb2^.  •  •  *  *  ^Sl^  * 

and  the  second  to  occur.  This  will  occur,  on  the  average,  earlier  with 
increasing  S  ,  so  the  approximations  are  even  better.  (11.1)  and  (11.2) 
are  first  modified  by  replacing  -  G^  by  the  sums  £  (H^  -  G^)  and 
h^  and  g^  become  h^  and  g^  ,  in  an  obvious  notation.  Then 
(11.1)  and  (11.2)  become: 

(11.1")  m(T)  a  (2/SyT2)  -  cxU/yt^)  , 

I (T)  a  S2y2T4/4  . 


(11.2") 


JU 

Not  only  is  at. .  <  (2/S)  the  new  requirement  that  u)(t...)  be  positive, 

*1  *1 

the  other  bound  becomes:  (2/ S) (1  -  (2y/oiS))  ,  thus  widening,  relatively, 
the  gap  in  which  t.,  must  fall  to  give  an  initial  estimate  corresponding 
to  reliability  growth!  The  quadratic  dependence  of  I  on  S  is  also 
reassuring. 

Unfortunately,  discussion  of  the  next  interval  of  T  ,  beyond  the 
second  failure,  already  leads  to  quadratic  expressions  for  w(T)  . 

C.  Average  Behavior  for  yT  Small  and  S  Large 

Before  giving  the  results  of  this  section,  let  us  recall  what  kind  of 
results  we  would  get  if  we  were  estimating  the  failure  rate  X  of  a 
stationary  Poisson  process,  with  S  systems  on  test  for  T  time  periods. 
The  maximum  likelihood  estimator  of  X  is  just  X (T)  =  N(T)/ST  ,  and  it 
is  immediate  that  this  estimator  is  unbiased  for  all  T  ,  since  E{ X (T) }  = 
SH (T) / ST  =  X  for  all  T  .  Furthermore,  E{ I (T) }  =  SM(T)  =  SXT  ,  so  that 
the  average  amount  of  precision  in  any  experiment  is  proportional  to  ST  , 
the  total-time-on- test . 

To  obtain  similar  results  using  (10.1),  using  h(t)  =  1  but  arbitrary 
g(t)  ,  we  first  of  all  set  T  small  enough  so  that  the  probability  of  more 


than  one  failure  per  system  in  (0,T]  is  negligible  compared  to  the  prob¬ 
ability  of  zero  or  one  failure.  If 


(11.4)  g(t)  *  1  -  y^t  +  y 2t2/2  “  Y3t3/6  +  ...  (t  -*•  0)  , 


Letting  ti»t2’  CN(t)  denote  c^e  random  number  of  (first) 

failures  which  would  be  observed  in  such  a  case,  we  rewrite  (10.1)  as: 

N(t)  1  -  g(t  ) 

(11.5)  T  -  G(T)  =  T  l  - - -  • 

k=l  ag (t^)  +  u(l  -  gCt^)) 

Now  the  t^  are  independent,  identically  distributed  random  variables 
with  the  same  distribution  as  the  {t  ^  |  t^  £  T}  ,  that  is,  they  have 
a  normalized  density  with  failure  rate  g(t)  that  is  approximated  by  (11.4) 
With  a  little  bit  of  labor,  we  find  the  first  two  moments  of  these  first 
failure  epochs  to  be: 


EUk}  =  |  - 


+  (a  -  ai)y 


12a 


-  T2  +  0(T3) 


(11.6) 


+  0(T3) 


Now  let  the  number  of  systems  on  test  get  very  large;  then  (11.5)  becomes: 


T  -  G(T)  =  E< 


1  -  g(tk) 


/ag(tk)  +  a) (1  -  g(tk))1 


almost  surely  as  S  -►  00  . 


We  now  expand  both  sides  of  this  equation  in  powers  of  T  ,  using  (11.4) 
and  (11.6),  assuming  that  E{a> (T) }  2;  +  0(T  )  .  After  some  labor, 

we  find  that  the  first  and  second  powers  of  T  are  trivially  satisfied, 
but  that  from  the  third  power  we  get: 

2 

E{w(T)}  ss  w  -  | —  +  0(T)  (YjT  «  1)  . 


(11.7) 


This  means  that  the  maximum  likelihood  estimator  is  inconsistent  at 


small  values  of  T — in  fact, for  the  parameters  of  the  numerical  trials, 
oi (T)  approaches  oj  -  5.0  !  This  gives  us  additional  evidence  as  to 
the  unsuitability  of  the  maximum  likelihood  estimator. 

Using  similar  techniques,  we  can  estimate  the  average  information 
in  the  data  for  small  T  as: 

2_3 
Y  T 

(11.8)  E(I(T)}ssS-^-+0{T4)  (YjT  «  1)  . 

Since  the  mean  number  of  samples  is  of  the  order  of  SaT  ,  this  shows 
that  the  precision  of  the  estimate  is  increasing  strongly  with  increasing 
time-on-test  T  ,  not  just  proportionally  to  total-time-on-test  ST  . 

For  the  data  in  the  numerical  trials,  this  estimate  of  precision 
is  ST  /1200  .  Table  I  shows  that  this  estimate  is,  in  fact,  rather  good 
for  yT  _<  1.0  ,  and  S  >_  4  ,  thus  explaining  the  regularity  observed  in 
Figure  9.  Remember,  however,  that  this  refers  to  the  second  derivative  of 
the  £n- likelihood  at  a  possibly  negative  i  . 


T 

T 

T 

T 


TABLE  I.  Data  Information,  I,  for  Different  S  and  T, 
Estimated  from  Numerical  Trials  (upper  numbers) 
and  from  Approximation  (11.8)  (lower  numbers). 


as 


D.  Bayesian  Estimators 


This  is  perhaps  a  convenient  point  to  re-emphasize  that  none  of  the 
above  difficulties  occur  when  using  a  Bayesian  analysis.  Not  only  does 
(9.2)  give  the  posterior-to-experiment  distribution  of  w  for  any  testing 
protocol,  it  follows  then  that,  given  a  risk  function,  we  can  make  point 
estimates  of  the  parameters,  or,  if  desired,  can  specify  Bayesian  tolerance 
regions  for  p(w  |  V)  .  If  either  the  prior  or  the  likelihood  is  approxi¬ 
mately  normal,  we  can  approximate  the  posterior-to-experiment  mean,  E(w  |  V ) 
by  a  precision-weighted  "credibility"  [12]  mixture  of  the  prior  mean,  E(uj)  , 
and  the  MLE  w  ;  the  precision  of  this  estimate  is  just  the  sum  of  the  prior 
precision  and  the  data  information.  And  so  on.  There  are  no  conceptual 
difficulties  with  estimating  more  parameters,  although,  as  indicated  above, 
numerical  integration  would  be  necessary  in  most  cases. 


12.  CONCLUSIONS 


Before  summarizing  the  results  of  this  study,  we  recall  the  basic 
assumptions  behind  our  learning-curve  reliability  growth  model :  only 
the  past  epochs  of  failure  may  be  used  to  estimate  process  parameters 
(although  the  modes  of  failure  may  be  used  in  making  the  actual  improve¬ 
ments)  ;  the  date  of  the  most  recent  failure  of  a  system  is  a  surrogate 
for  all  previous  history  of  failures  and  modifications  of  that  system, 
that  is,  the  system  state  is  Markovian  at  these  epochs  and  depends  only 
upon  clock  time  and  local  age;  and,  the  age  and  global  time  failure  rate 
effects  are  separable,  as  in  (2.2). 

Given  this  framework,  we  believe  the  following  conclusions  can  be 

made: 

(1)  A  complete  learning  curve  reliability  growth  model  requires 
at  least  two  parameters  and  two  prototypical  failure  rate 
functions  to  describe  the  desired  phenomenon; 

(2)  If  the  learning-curve  effect  is  continuous  and  slow,  the  exact 
form  of  the  prototypical  growth  function  g  is  probably  not 
too  important;  there  is  also  not  much  difference  between  the 
"as  operated"  and  "as  produced"  specializations  introduced 

in  Section  3; 

y— 1 

(3)  However,  use  of  the  Duane  learning  curve,  g(t)  =  Kt  ,  leads 

to  technical  difficulties  in  the  problems  of  interest.  An  ex- 
ponential  learning  curve,  g(t)  ■  e  ,  avoids  these  difficulties, 
and  is,  moreover,  suggested  by  certain  defect  removal  models  of 
software  reliability  growth; 

(A)  In  any  case,  the  mathematical  problems  of  finding  the  maximum- 
likelihood  estimates  of  the  parameters  in  (6.1),  (6.2),  (10.1), 


or  the  data  information  (matrix)  measure  (6.3),  (6.4),  (6.5), 
(10.2),  are  similar  for  the  different  models  and  forms; 

(5)  Limited  computational  experience  in  estimating  the  ultimate 
failure  rate,  to  ,  under  favourable  assumptions  shows  that 
the  maximum-likelihood  estimator,  ui  ,  is  very  unstable 
during  the  early  part  of  the  learning  curve  (y^T  <  1)  , 
especially  for  a  small  number  S  of  systems  on  test.  In 
fact,  for  y^T  small  and  S  large,  to  is  inconsistent, 
usually  much  less  than  the  true  to  ,  and  negative; 

(6)  The  estimate  of  the  data  information  measure,  on  the  other 
hand,  is  remarkably  stable.  (11.8)  is  apparently  valid  as 
a  scaling  law  for  y^T  1  and  S  >_  4  ,  and  can  be  used  to 
estimate  the  precision  of  relatively  non- informative  testing 
protocols.  In  any  case,  it  is  always  better  to  increase  T 
than  to  increase  S  in  this  region. 

Thus,  the  estimation  of  ultimate  performance  using  classical  estima¬ 
tion  methods  presents  fundamentally  difficult  problems,  unless  testing 
periods  are  extended  well  into  the  learning  curve,  and  a  reasonable  number 
of  systems  are  placed  on  test.  The  problem  of  simultaneous  estimation  of 
several  parameters  will  likely  introduce  additional  inconsistencies  and 
instabilities. 

Bayesian  estimation  procedures,  on  the  other  hand,  use  data  from  any 
testing  protocol  in  a  consistent  and  "friendly"  manner,  with  the  posterior- 
to-testing  density  of  (o  given  directly  from  the  likelihood  and  the  prior. 
Point  estimates  or  Bayesian  tolerance  regions  can  then  be  obtained  through 
numerical  integration,  or  through  approximations,  if  either  the  prior  or 
the  likelihood  has  high  precision.  We  expect  to  see  further  development 
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